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Rejection-free Monte Carlo Algorithms for Models with Continuous Degrees of 

Freedom 
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We construct a rejection-free Monte Carlo algorithm for a system with continuous degrees of 
freedom. We illustrate the algorithm by applying it to the classical three-dimensional Heisenberg 
model with canonical Metropolis dynamics. We obtain the lifetime of the metastable state following 
a reversal of the external magnetic field. Our rejection-free algorithm obtains results in agreement 
with a direct implementation of the Metropolis dynamic and requires orders of magnitude less 
computational time at low temperatures. The treatment is general and can be extended to other 
dynamics and other systems with continuous degrees of freedom. 
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Nucleation and metastability are characteristic behav- 
iors of dynamical processes for many different fields, from 
stock markets and sociology Q to parallelization meth- 
ods for massively parallel computers || to chemical reac- 
tions and materials science. Many classical models, when 
simulated with Monte Carlo methods, also present these 
behaviors, and although the number of steps in the sim- 
ulation does not necessarily correspond directly to ex- 
perimental time, they give valuable insights into these 
dynamic phenomena. For instance, studies on Ising ||] 
and anisotropic Heisenberg models Q have shown the 
existence of different metastable decay regimes for small 
ferromagnetic particles after a reversal of the external 
magnetic field. At low temperatures, single and multi- 
ple droplet nucleation and a strong field regime are ob- 
served, and recently, indirect experimental evidence of 
these regimes has been found [|| . 

Many Monte Carlo dynamics are Markov processes 
that divide each step into two successive parts: first, a 
new state is chosen, second, it is accepted or rejected 
according to some criteria. In many cases of interest, 
the acceptance rates can be so small that a huge num- 
ber of trials is required before a change is made. Then, 
a direct implementation of the Monte Carlo dynamic, 
one that attempts steps one after the other, is extremely 
slow. For instance, Ising models with Metropolis dynam- 
ics H can require 10 15 trials to leave a metastable state 
at low temperatures, and such a simulation would take 
10 10 minutes &. Therefore, techniques to implement the 
same dynamic in a faster way are required. 

There exist different techniques to construct rejection- 
free implementations of Monte Carlo dynamics for dis- 
crete spin systems, like the n-fold way Jg|, ^ and Monte 
Carlo with absorbing Markov chains [^[algorithms (for 
a review see fllO|). In this paper, we extend the n-fold 



way rejection- free technique to systems with continuous 
degrees of freedom, and we construct a rejection-free al- 
gorithm for the classical Heisenberg model. Our treat- 
ment is completely general and can be extended to other 
dynamics and other continuous systems. 

Consider a Markov process with every step consisting 
of two parts. First, choose a movement from state x to 
state x' ^ x with probability T(x'\x). Second, accept it 
with probability A(x'\x). The full probability to undergo 
the movement {x'\x) is then 



W{x'\x) = T(x'\x)A{x'\x) 



(1) 



In this paper, the term "direct implementation" refers 
to the common selection and rejection implementation of 
the Markov process, in which two random numbers are 
used, one for the selection of x' and one for the rejec- 
tion/acceptance of x' . 

In the specific case of importance sampling, T(x'\x) is 
chosen to be symmetric, T(x'\x) = T(x\x'), and A{x'\x) 
is tuned to obtain a desired limit probability distribution, 
P(x), when the number of steps tends to infinity. This is 
accomplished by requiring the detailed balance condition, 



W(x'\x)P(x) = W(x\x')P(x') 



(2) 



To obtain the canonical distribution, P{x) cx 
exp[— E x /ksT], a widely used choice for A(x'\x) is the 
Metropolis acceptance probability 



A(x'\x) = vmin{l,exp[-(E x , - E x )/k B T}} 



(3) 



where v is a constant which is the same for all x and 
x' and is included to ensure that J2 X ' W(x'\x) < 1 for 
all algorithmic steps. For all discussions presented here, 
v = 1 was sufficient to ensure this condition. The usual 
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Metropolis simulation g is a direct implementation of 
this dynamic. 

The n-fold way |§, 0] is a rejection-free implementa- 
tion of the dynamic W(x'|x), and we briefly remind the 
reader of this two-part algorithm. First, the number 
of trials, t, to leave the current state is computed (the 
update time), and second, one movement is chosen and 
performed. In this way, every algorithmic step induces 
a change in the system configuration, but the dynamic 
information is preserved. 

The first step of the n-fold way algorithm consists of 
generating the update time, t, which is a random variable 
chosen from the appropriate probability distribution. De- 
fine A as the probability to reject all movements, 



X = 1-J2 W ( X '\ X ) > 



(4) 



and the probability p(t) to leave state x after t steps is a 
geometric distribution 0, 



o(t)=X t - 1 (l-X) 



(5) 



A so-called integral generator [[fl| can be constructed 
to produce t with this distribution. Define I(t):=l — 
^ fe=1 p(fc)=A* (I(0):=1), and let f be a random number 
uniform on (0,1). The number of steps until the next 
update, t, is then determined by 



I t -x < r < I t and t 



lnf 
m~A 



(6) 



where [x\ denotes the integer part of x. 

In the second step of the n-fold way algorithm, we must 
choose the exit state, x' , from the appropriate probability 
distribution. The probability to exit to state x' from x 
is 



C{x'\x) 



1 



1-A 



W{x'\x) 



(7) 



A movement is chosen with this probability by using 
the same integral-generator strategy as in the first step 
of the algorithm. The movements are ordered with in- 
dex k (discrete degrees of freedom), and the partial 

sums Q(k):=^2 m —i C(x'\x), n , with Q(0):=0, arc com- 
puted. The movement k is chosen if Q(k — l)<f<Q(k), 
where f is uniformly distributed on (0,1). To per- 
form the movement selection in a more efficient way, 
movements are grouped into classes, and a two-level 
search is performed by first selecting the class i with 
probability C(i\x)=^2 x , ei C(x'\x) and then selecting the 
movement (x'\i) from inside the class with probability 
C(x'\i)=C(x'\x)/C(i\x) (Bayes). This completes the n- 
fold algorithm, which is a rejection-free Monte Carlo al- 
gorithm for systems with discrete degrees of freedom. 

In importance sampling for discrete spin systems J^, H . 
the exit states, x', are usually grouped into classes by 
energy changes, where all movements in class i have the 
same energy change, AE, and thus the same acceptance 



rate, A4. Let be the number of movements in class 
i. Using T(x'\x)=l/N, where N is the total number 
of possible movements, we obtain the classic n-fold way 
expressions M 



X = 1 _ y V± Ai ; C(i 

; N 



TliAi 

(1 - X)N 



c(x'\i) 



1 

n, 



(8) 



To apply this rejection-free technique to systems with 
continuous degrees of freedom, some ideas employed to 
extend the broad histogram method for continuous sys- 
tems are used @ . Discrete probability distributions be- 
come probability density functions (pdf), and a discrete 
choice of a random variable becomes the construction 
of a random generator with its proper distribution. To 
illustrate this process, we construct a rejection-free algo- 
rithm for the classical Heisenberg model with Metropolis 
dynamics. The Hamiltonian is 



H 



{JxXjXj + JyYiYj + J z ZiZj} — H z 



2>. 



where Ui=XiX + Yiy + Ziz is a spin of unit length on site 
i, H z is the magnitude of an external magnetic field in 
the z-direction, (ij) represents a nearest- neighbor sum- 
mation, and J Xl J y , and J z are coupling constants. 

For continuous systems like the Heisenberg model, the 
movements form an uncountable set. With x fixed, 
T(x'\x) can be interpreted as a pdf in the configuration 
space of the system, where T(x'\x)dx' is the probabil- 
ity to choose the new state, x', inside a small infinites- 
imal region dx' , centered on x' . Let us fix T(x'\x) by 
choosing all movements with the same probability as fol- 
lows. First, a site i is chosen with probability l/N, where 
N is the number of sites. Second, a new orientation a[ 
for the spin at i is chosen uniformly on the unit sphere 
(pdf T(9',ip'\i)—l/4Tv). This is equivalent to generating 
z'=cos9' uniformly on the interval [—1,1] and ip' uni- 
formly on the interval [0, 2tt), where (z 1 , if') are the coor- 
dinates of a[ in some cylindrical coordinate system [[12] . 
Let (z, ip) be the coordinates of <7j in the same system . 
The total pdf of this movement (z' , <p'\z, ip)i is, therefore, 
T(z',(f'\z,(p)i=l/4TrN, and clearly, T{x\x')=T{x'\x). 

The energy change of this movement is 



(10) 



Si = 





AE-- 


= & - ?■) 


■Si , 


Jx Xj 


x + 




y + 
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Here, j denotes a sum over the nearest neighbors of site 
i. If we rotate to a coordinate system with the z-axis 
parallel to Si, this energy change reduces to AE=—(z' — 
z)S l . Therefore, from Eq. (3), 



A(z',ip'\z,ip), 



_ J exp[(z' 



i)Si/k B T] 



if z < z 
otherwise 



(11) 
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Thus, W(z', <p'\z, v )i=(l/4nrN)A(z', <p?\z, tp)t. 

To implement the rejection-free algorithm, we start by 
computing A for this dynamic 



z' such that Qi(z')- 



1 N 



(12) 



i=l 



Xi := f* = _ 1 f^ =0 T(z',<p'\i)[l-A(z',(p'\z,<p) i ]dip'dz' 
= ^ - W t 1 _ ex P(-(^ + W*bT)] . 

Finally, the value of i to exit from the current state is 
computed from Eq. (6). 
According to Eq. (7), 



C{z',<p'\z,<p)i 



4ttN{1 - A) 



A(z',(p'\z,<p)i 



(13) 



but for a system with continuous degrees of freedom, it 
is not possible to group the movements into classes by 
energy changes, since these values form a continuous set, 
and instead of grouping by energies, we group the move- 
ments by sites. The probability to choose a site i is 



C{i\x) 



C(z',ip'\z,ip)i dz'dtp' = J. \ 

sphere JV K 1 ~ A ) 



The search to find which spin i changes can be per- 
formed with the same integral-generator strategy de- 
scribed above. 

Next, we choose (z 1 , (p 1 ) for the site i, but since these 
are continuous variables, their distribution is described 
by a pdf that, according to the Bayes relation, is given 
by 



C{z',ip'\i) = 



1 



4tt(1 - Xi) 



A(z\<p'\z,<p)i 



(15) 



Since this expression is independent of ip', this coordi- 
nate is uniformly distributed on the interval [0, 2n). In 
contrast, z' must be generated with pdf 



P.-'; - J 2(1-A S ) 
2(1-Ai) 



exp((z' - z)S % /k B T) ; if z' < z 
: otherwise 



(16) 

A rejection-free random generator with this distribu- 
tion can be constructed by means of the same integral- 
generator strategy but on a continuous variable. First, 
define the integral 

QiCO := j:L- 1 fi(*")dz" 

2 { i-f i)Si ( ex P((^' - z)Si/k B T) - 1] + tU(z) 
; if z' < z 



otherwise 



(17) 

with Qi(z)=l — (1 — z)/(2(l — Xi)). Next, generate a ran- 
dom number f uniformly distributed on (0, 1) and take 



a + ^-ln [2(1 - Xi)Si(f 
z + 2(1 — Ai)(f — £li(z)) ; otherwise 



ni(z))/k B T+l] 
if f < Qi{z) 



Finally, the new values {z' , ip') are transformed back into 
Cartesian coordinates and rotated back to the original 
reference axis. Thus both the update time, t, and the exit 
state, x' , are chosen, and the rejection-free implementa- 
tion for a system with continuous degrees of freedom is 
complete. 

We tested our rejection- free algorithm by computing 
the mean lifetimes of a metastable state for an anisotropic 
Heisenberg model with J x =J y =l.O and J z =2.0 on a 
simple cubic (SC) lattice of size 10x10x10 with peri- 
odic boundary conditions. The system is initially in a 
metastable state with all spins pointing in the —z direc- 
tion with the external field in the opposite direction, and 
all simulations were performed at temperatures T < T c , 
where we have found the approximate value of the critical 
temperature to be T C ~3.15J X . Here, we concentrate on 
dynamic quantities, such as the metastable escape time 
or lifetime, r, which is the number of Monte Carlo steps 
per site (MCSS) needed to obtain a magnetization M z =0. 
Preliminary results from the algorithm for static quanti- 
ties have been presented elsewhere [ p^| . 

To improve the efficiency of the site selection portion of 
the algorithm, the sites were grouped into super-classes 
of lines and planes and organized into a three-level tree 
Q. Even with this improved tree-search, the steps of 
the rejection-free algorithm require more computational 
time that the steps of a direct Metropolis implementa- 
tion. On average, the rejection-free method took 8.769/is 
to perform one change, and the direct Metropolis im- 
plementation took 1.162/is to make one trial. However, 
when many trials are required to make a single change, 
the rejection- free algorithm is computationally more ef- 
ficient than the direct implementation. 

Figure [l] shows the average lifetimes, (r), of the 
metastable state computed by both rejection-free and 
direct implementations vs. 1/T for many values of the 
external magnetic field. As expected, the two methods 
give identical results for all temperatures and fields, re- 
gardless of the switching mechanism (nucleation regime 
H z <6.0, strong-field regime H z >6.0). See Ref. jl5| for a 
description of the switching mechanisms. 

We define (t)cpu, direct as the average CPU time re- 
quired to escape the metastable state, when simulated 
by the direct implementation, and (T}cpu,rej-free is sim- 
ilarly defined. Figure |^ shows the CPU-time ratio, 

(T)cPU,direct/(T)cPU,rej-free, VS. 1/T. For T = 1/100, 

the rejection- free implementation is nearly two orders of 
magnitude faster than the direct implementation. 

In summary, we have constructed a rejection- free 
Monte Carlo algorithm for a system with continuous de- 
grees of freedom which faithfully keeps the dynamic as 
compared to the direct implementation but is orders of 
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magnitude faster at low temperatures. Our procedure is 
quite general and can be applied to other dynamics and 
other systems with continuous degrees of freedom. 



1 10 100 



1/T 



FIG. 1: Average lifetimes, (r), of the metastable state of 
an anisotropic Heisenberg model on a 10 3 cubic lattice at 
many magnetic field values. Each point represents an av- 
erage over fOO independent escapes. Results for both the 
rejection-free algorithm (triangles) and the direct implemen- 
tation (filled squares) are shown. The vertical line represents 
the approximate value of T c . 

100 r ■ ■ ■ I 



1 10 
CD 

E 
7 

b 1 



1 10 100 

1/T 



FIG. 2: CPU time ratio vs. 1 /T. The external field values are 
H z —5.6 (squares), H z —5.96 (circles), and H z =7.0 (triangles). 
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